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Abstract. The problem of equilibrium of a plasma in a Tokamak is a free boundary problem 
described by the Grad-Shafranov equation in axisymmetric configurations. The right hand side 
of this equation is a non linear source, which represents the toroidal component of the plasma 
current density. This paper deals with the real time identification of this non linear source from 
experimental measurements. The proposed method is based on a fixed point algorithm, a finite 
element resolution, a reduced basis method and a least-square optimization formulation. 



1. Introduction 

In Tokamaks, a magnetic field is used to confine a plasma in a toroidal vacuum vessel. The 
magnetic field is produced by external coils surrounding the vacuum vessel and a current 
circulating in the plasma, forming a helicoidal resulting magnetic field. At equilibrium, magnetic 
field lines lie on isosurfaces forming a family of nested tori and called magnetic surfaces. These 
magnetic surfaces enable to define the magnetic axis and the plasma boundary. The innermost 
magnetic surface which degenerates into a closed curve is called magnetic axis. The plasma 
boundary corresponds to the surface in contact with a limiter or being a magnetic separatrix 
(hyperbolic line with an X-point). 

Denote by j the current density, B the magnetic field and p the kinetic pressure. The equilibrium 
of the plasma in presence of a magnetic field is described by 

j = Vx-, (1) 

jxB = Vp, (2) 
V-B = 0, (3) 

where /x is the magnetic permeability. Equation ([T|) is Ampere's theorem and equation ([3]) 
represents the conservation of magnetic induction. Equation ([2]) means that at equilibrium the 
Lorentz force j x B balances the force Vp due to kinetic pressure. From this equation it is clear 
that 

B • Vp = and j • Vp = 0. (4) 



Thus, p is constant along magnetic field lines and current lines which consequently lie on 
magnetic surfaces. 



Figure 1. Toroidal geometry. 



Consider the cylindrical coordinate system {r,z,(p). Under axisymmetrical hypothesis, the 
magnetic field is supposed to be independent of the toroidal angle 4> and can be decomposed in 
the form 

B = Bp + B J' 

Bt = -60 

where Bp and By are respectively the poloidal and the toroidal component of the magnetic 
field, e^f) is a unit vector in direction (j),the term 'i/;(r, z) is the poloidal magnetic flux function 
(see Figure [1]) and / is the poloidal current flux function. The magnetic surfaces are defined 
by rotation of the flux lines tp = constant around the vertical axis (Oz). Using ([5]), equations 
dl])-® lead to the Grad-Shafranov equation (see [HEIIS]) 

_^*^ = ^p'^^^ + J_^ff'^^^^ (6) 

where 

did did 
A*. = —(——) + —(——). (7) 

dr fior dr dz fi^r dz 

and ^0 is the magnetic permeability of the vacuum. Thus, under axisymmetric hypothesis, the 
three dimensional equilibrium ([l])-([3]) reduces to solve a two dimensional non- linear problem. 
Note that the right hand side of ([6]) represents the toroidal component of the plasma current 
density which is determinated by p', / and /'. 

In this paper, we are interested in the numerical reconstruction of the plasma current density 
and of the equilibrium ([6]) (see [H El [6]). This reconstruction has to be achieved in real time 
from experimental measurements. The main difficulty consists in identifying the functions p' 
and //' in the non linear source term of ([6]). An iterative strategy involving a finite element 
method to solve the direct problem ^ and a least square optimisation procedure to identify the 
non linearity using reduced basis is proposed. A description of the experimental measurements 
available in Tokamaks is given in Section 2. Section 3 is devoted to the statement of the 
mathematical problem and the numerical algorithm proposed. Numerical results obtained with 
the software Equinox are presented in Section 4. 



2. Experimental measurements 

Although p and / cannot be directly measured in a Tokamak, several measurements are available 
(see Figure [2]) : 



magnetic measurements: the flux loops provide Dirichlet condition ip = h and magnetic 

1 dip 

probes Neumann condition — — — = g on the walls of the vacuum vessel which represent the 

r on 

boundary of the computational domain. These measurements are given at some discrete 
points Ni and the Dirichlet full boundary data are obtained by interpolation. 

polarimetric measurements: this measurement gives the value of the integral along a family 
of chords Cj 

r ne dip ^ 
Jci r dn 

dip 

neiip) is the electronic density which is approximately constant on each flux line, is the 

dn 

normal derivative of ip along the chord Q. 

interferometric measurements: they give the values of the integrals / n(.dl = (5i. 
current measurements: they give the value of the total plasma current /„ deflned by 



-dl = ai. 



j^dx. 



Other measurements are potentially available but are not used for the moment in the software 
Equinox. 




Figure 2. Left: the straight green lines represent the chords used for polarimetry and 
interferometry measurements. Right: part of the vacuum vessel including the magnetic 
measurements. At the bottom middle an example of flnite element mesh used for numerical 
simulations. 



3. Algorithm and numerical resolution 

Let Q be the domain representing the vacuum vessel of the Tokamak, and dQ its boundary. The 
equilibrium of a plasma in a Tokamak is a free boundary problem. The plasma boundary is 



determinated either by its contact with a limiter D or as being a magnetic separatrix with an 
X-point (hyperbohc point). The region Qp C containing the plasma is defined by 

$7p = {x G U, '0(x) > tpb} 

where iph = max£) ip in the hmiter configuration or ipi, = il>{X) when an X-point exists. 
In the vacuum region, the right hand side of ([6|) vanishes and the equihbrium reads 

A*ip = in O \ Op 

ijj — max ip 

Consider the fohowing notations: = € [0, ll in Qp , j4(^) = ——p'ii^) and 

ipb — max ip X 

-^(V') — ~{ The functions A and B are fiux functions defined on the fixed interval 

[0, 1]. Giving Dirichlet boundary conditions, the final equilibrium problem is 

Ro r " (8) 

■0 = /i on do, 

where is the characteristic function of 0,p and Rq is the major radius of the Tokamak. The 
parameter A is a normalization factor that satisfies 

Ip = x[ [^Aii,) + ^B{^|J)]dn. (9) 

3.1. Iterative algorithm 

The aim of the method is to reconstruct the plasma current density and the equilibrium solution 
in real time. At each time step determined by the availability of new measurements, the 
method consists in constructing a sequence (V'", ^^p, -B") converging to the solution vector 
(V', f^p, A, B). The sequence is obtained by the following algorithm: 

• Starting guess: ip^ , Op, AP and B^ known from the previous time step solution. Compute 
A satisfying ([9]). 

• Step 1 - Optimisation step: computation of using a least square 
procedure. 

• Step 2 - Direct problem step: computation of t/;""''^ and Op"*"^ solution of 

-A*0-+i = X[^A^^\r) + ?^B^+\r)]xn-, inO 

^0 (10) 

0"+i = /i on ao. 

• Step 3: If the process has not converged, n := n + 1 and return to Step 1 else 
{tpjAjB) = ('0"', yl", S"). The process is supposed to have converged when the relative 

residu rr- — n is small enough. 

At each iteration of the algorithm, an inverse problem corresponding to the optimization step 
and an approximated direct Grad-Shafranov problem corresponding to pi?|) have to be solved 
successively. In (jlOp . ip^ is supposed known, the right hand side does not depend on ip^^^ and 



thus problem (jlOp is linear. 

Step 1 consists in a least-square minimization formulation 

Find A*, B*, n* such that : 

(11) 

J{A*,B*,nl) = inf J{A, B, n^). 

When polarimetric and interferometric measurements are used, the electronic density has also to 
be identified even if rig, depending on ^, does not appear in ([5]). The cost function J is defined 
by 

J(A,S,ne) = Jo + i^lJl +i^2J2 + (12) 

where 

^0 = Y.^-^-lm - h = E( / ^^dl - a.)^ and J, = E( / n^dl - A)^ 

i i i 

and the coefficients Ki and K2 are weights giving more or less importance to measurements used 
(see [5j). As a consequence of the ill-posedness of the identification of A, B and ng, a Tikhonov 
regularization term is introduced (see [7]) where 

J, = ei (^[A"{x)fdx + e2 [B" {x)f dx + (^K{x)fdx 
Jo Jo Jo 

and ei, £2 and are regularization parameters. 

In the next section, the numerical methods used to solve the two involved problems are detailed. 



3.2. Numerical resolution 

The resolution of the direct problem is based on a finite element method (see [8j). Consider 
the family of triangulation of Q, and the finite dimensional subspace of H^{^1) defined by 

Vh = {vh G H\Q),Vh\T G P\T), Vr G rj. 

Introduce = Vh H Hq{Q), the discrete variational formulation of PU]) reads 

Find iph €z Vh with iph = h on dQ such that 

' n f I f r - Rn - (13) 

yvh G V^, / • Vvhdx = / X[—A{r) + —B{r)]vhdx 

Jn nor jQ,p Ro r 

where ■0* represents the value of -0 at the previous iteration. The Dirichlet boundary conditions 
are imposed using the method consisting in computing the stiffness matrix of the Neumann 
problem and modifying it. The modifications consist in replacing the rows corresponding to 
each boundary node setting 1 on the diagonal terms and elsewhere. The Dirichlet conditions 
appear in the right hand side of the linear system. At each iteration, only the right hand side 
of the system has to be modified. The sparse matrix K and its inverse are computed at 
the beginning of the algorithm and stored until the end of the simulation. 
The linear system of (|13p can be written in the form 

K."^ = y + H (14) 



where K is the modified stiffness matrix, ^' is the unknown vector, y is the right hand side of 
the problem and H is the term corresponding to the Dirichlet conditions. The vector y depends 



on A and B determinated in the optimization step. 

The functions A, B and Ue are decomposed on a basis (<5i)j=i, °f dimension m 

mm m 

A{x) = ai^i{x), B{x) = hi^i{x) and ne{x) = ^ Ci^i{x). 

i i i 

The vector y reads 

y = (15) 

where u = (oi, am, ^m) £ R^™. The term y is a matrix of size n x 2m where n is the 
number of nodes. Consider (f j) a basis of Vh, each row z of y is decomposed as 

/ \—^;-rn{V)vd^ if m + 1 < j < 2m. 

During the optimisation step, rig is first estimated from interferometric measurements and A 
and B are computed in a second time. The function neiijj) is approximated using a least square 
formulation for the minimun of J2 with Tikhonov regularization, and solving the associated 
normal equation. 

To approximate A and B, suppose He is known and consider the discrete approximated inverse 
problem 

Find u minimizing : 

J{u) = \\\C{r,ne)^ - k\\l + |n^Au ^^^^ 

where C{^*,ne) is the observation operator corresponding to experimental measurements given 
in Section 2 and k represents the experimental measurements. The first term in J is the discrete 
version of Jq + i^i Ji- The second one corresponds to the first two terms of the Tikhonov 
regularization term J^. Denote by / the number of measurements available and by af the 
variance of the error associated to the i^^ measurement, the norm \\.\\d is defined by 

Vx G Ml = (Z)x,x) = (1)1/2^, Z?i/2x) 
where D is the diagonal matrix da = 

C is a matrix of size / x n and can be viewed as a vector composed of two blocks corresponding 
respectively to Jq and Ji. 

The matrix A is of size 2m x 2m and is block diagonal composed of two blocks Ai and A2 of 
size m X m, with 

(Ai)i, = (A2)y = / $",(x)$",(x)dx 
Jo 

where the $j are the second derivatives of the basis functions <I>j. 
Using (fT^ and (fTSj) . the problem (fT6]) becomes 

J(u) = ]^\\C{^*,ne)^ -k\\% + UjL^Ku 

= ]^\\C{r,ne)K-^Y{r)u + {C{r,n,)K-^H - k)\\l + ku 



= ^\\Eu-F\\l + ^u^Au 

where E = C{i>* ,ne)K-'^Y{ip*) and F = -C{ij* ,ne)K~'^ H + k. Setting E = D^/^E, problem 
(|16|) reduces to solve the normal equation 

{E^ E + eA)u = E'^ F (17) 



Table 1. Numerical convergence. 
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7.43857e-07 
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4.1659e-06 



4. Numerical results 

The method detailed in this paper has been implemented in the software Equinox developed 
in collaboration with the Fusion Department at Cadarache for Tore Supra (see [9} 110^ [TT]) and 
JET (Join European Torus). Equinox can be used on the one hand for precise studies in which 
the computing time is not a limiting factor and on the other hand in a real-time framework. 

In table [H the evolution of the relative residu onip^ A and B versus the number of iterations 
is given. It demonstrates numerically the convergence of the algorithm. 

If a value e = 10~^ is used as stop condition, the algorithm almost always needs about 5 to 
30 iterations to converge. For the finite element mesh of 412 nodes and 762 elements used at 
Tore Supra the average iteration cost (file access included) is of 0.038 s on a processor Intel(R) 
Core(Tm)2 Duo T7100 l.SOGHz. The expensive operations are the updates of matrices C and 
Y and the computation of products CK~^ and CK~^Y (see Section [3. 2p . The resolution of the 
direct problem (Eq. [T^ is cheap since the matrix is stored once for all and does not change 
from one iteration to the other. Morevover the resolution of the normal equation (Eq. [T7j) is also 
cheap. Each function A and B is decomposed in a basis (which can be chosen to be cubic 
B-splines, piecewise linear functions or wavelet scaling functions) of about 5 to 10 functions 
(depending on the user's choice). This makes the dimension of the matrix to be inverted at 
most 20 X 20. 

For real-time applications the numerical reconstruction of an equilibrium must not take more 
than 0.050 to 0.100 s. Therefore in this context it is not possible to let the algorithm fully 
converge. A maximum number of iterations has to be set (typically 1 or 2 depending on the 
computer and the size of the mesh). In practice this is not a real problem since our approach 
is quasi-static. Two successive equilibriums during a pulse are very close and no important 
differences have been observed between the results of the fully converged algorithm and its 
real-time version. 

Concerning Tikhonov regularization with a real-time constraint it is absolutely not possible 
to compute the regularization parameter for each equilibrium. Therefore it is kept constant 
during a pulse. The L-curve method [12] provided a value of approximately 5 x 10~^. 

An example of the outputs of Equinox is presented in figure [3l It is a Tore Supra pulse in 
which magnetic, interferometric and polarimetric measurements on 5 chords are used. One can 
observe the position of the plasma in the vacuum vessel. Isofiux lines are displayed from the 
magnetic axis to the boundary. For each active chord, the error between computed and measured 



interferometry is given in purple. These errors are less than 1%. The polarimetry errors are given 
in yellow. They are one order of magnitude larger than for interferometry. Different graphs are 
plotted on the left of the display. On the first row the identified function A, and corresponding 
functions p' and p. On the second row the identified function B and corresponding function //'. 
The third row gives the toroidal current density in the equatorial plane and the fourth one 
shows the safety factor. Finally on the fifth row the identified rif, function is plotted. 




Figure 3. Example of Equinox output. See text for details. 
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